home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / r250.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  4.5 KB  |  169 lines
  1. /* rng/r250.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* This is a shift-register random number generator. The sequence is
  25.  
  26.    x_n = x_{n-103} ^ x_{n-250}        ("^" means XOR)
  27.  
  28.    defined on 32-bit words.
  29.  
  30.    The first 250 elements x_1 .. x_250 are first initialized as x_n =
  31.    s_n, where s_n = (69069*s_{n-1}) mod 2^32 and s_0=s is the
  32.    user-supplied seed. To ensure that the sequence does not lie on a
  33.    subspace we force 32 of the entries to be linearly independent.  We
  34.    take the 32 elements x[3], x[10], x[17], x[24], ..., 213 and apply
  35.    the following operations,
  36.  
  37.    x[3]   &= 11111111111111111111111111111111
  38.    x[3]   |= 10000000000000000000000000000000 
  39.    x[10]  &= 01111111111111111111111111111111
  40.    x[10]  |= 01000000000000000000000000000000 
  41.    x[17]  &= 00111111111111111111111111111111
  42.    x[17]  |= 00100000000000000000000000000000 
  43.    ....      ...
  44.    x[206] &= 00000000000000000000000000000111
  45.    x[206] |= 00000000000000000000000000000100 
  46.    x[213] &= 00000000000000000000000000000011
  47.    x[213] |= 00000000000000000000000000000010 
  48.    x[220] &= 00000000000000000000000000000001
  49.    x[220] |= 00000000000000000000000000000001 
  50.  
  51.    i.e. if we consider the bits of the 32 elements as forming a 32x32
  52.    array then we are setting the diagonal bits of the array to one and
  53.    masking the lower triangle below the diagonal to zero.
  54.  
  55.    With this initialization procedure the theoretical value of
  56.    x_{10001} is 1100653588 for s = 1 (Actually I got this by running
  57.    the original code). The subscript 10001 means (1) seed the
  58.    generator with s = 1 and then do 10000 actual iterations.
  59.  
  60.    The period of this generator is about 2^250.
  61.  
  62.    The algorithm works for any number of bits. It is implemented here
  63.    for 32 bits.
  64.  
  65.    From: S. Kirkpatrick and E. Stoll, "A very fast shift-register
  66.    sequence random number generator", Journal of Computational Physics,
  67.    40, 517-526 (1981). */
  68.  
  69. static inline unsigned long int r250_get (void *vstate);
  70. static double r250_get_double (void *vstate);
  71. static void r250_set (void *state, unsigned long int s);
  72.  
  73. typedef struct
  74.   {
  75.     int i;
  76.     unsigned long x[250];
  77.   }
  78. r250_state_t;
  79.  
  80. static inline unsigned long int
  81. r250_get (void *vstate)
  82. {
  83.   r250_state_t *state = (r250_state_t *) vstate;
  84.   unsigned long int k;
  85.   int j;
  86.  
  87.   int i = state->i;
  88.  
  89.   if (i >= 147)
  90.     {
  91.       j = i - 147;
  92.     }
  93.   else
  94.     {
  95.       j = i + 103;
  96.     }
  97.  
  98.   k = state->x[i] ^ state->x[j];
  99.   state->x[i] = k;
  100.  
  101.   if (i >= 249)
  102.     {
  103.       state->i = 0;
  104.     }
  105.   else
  106.     {
  107.       state->i = i + 1;
  108.     }
  109.  
  110.   return k;
  111. }
  112.  
  113. static double 
  114. r250_get_double (void *vstate)
  115. {
  116.   return r250_get (vstate) /  4294967296.0 ;
  117. }
  118.  
  119. static void
  120. r250_set (void *vstate, unsigned long int s)
  121. {
  122.   r250_state_t *state = (r250_state_t *) vstate;
  123.  
  124.   int i;
  125.  
  126.   if (s == 0)
  127.     s = 1;    /* default seed is 1 */
  128.  
  129.   state->i = 0;
  130.  
  131. #define LCG(n) ((69069 * n) & 0xffffffffUL)
  132.  
  133.   for (i = 0; i < 250; i++)    /* Fill the buffer  */
  134.     {
  135.       s = LCG (s);
  136.       state->x[i] = s;
  137.     }
  138.  
  139.   {
  140.     /* Masks for turning on the diagonal bit and turning off the
  141.        leftmost bits */
  142.  
  143.     unsigned long int msb = 0x80000000UL;
  144.     unsigned long int mask = 0xffffffffUL;
  145.  
  146.     for (i = 0; i < 32; i++)
  147.       {
  148.     int k = 7 * i + 3;    /* Select a word to operate on        */
  149.     state->x[k] &= mask;    /* Turn off bits left of the diagonal */
  150.     state->x[k] |= msb;    /* Turn on the diagonal bit           */
  151.     mask >>= 1;
  152.     msb >>= 1;
  153.       }
  154.   }
  155.  
  156.   return;
  157. }
  158.  
  159. static const gsl_rng_type r250_type =
  160. {"r250",            /* name */
  161.  0xffffffffUL,            /* RAND_MAX */
  162.  0,                    /* RAND_MIN */
  163.  sizeof (r250_state_t),
  164.  &r250_set,
  165.  &r250_get,
  166.  &r250_get_double};
  167.  
  168. const gsl_rng_type *gsl_rng_r250 = &r250_type;
  169.